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1.  Introduction 

Due  to  their  simplicity,  low-thrust  cold-gas 
propulsion  systems  have  the  potential  to  provide  or¬ 
bital  maneuvering  capability  for  small,- micro-  and 
even  nano-satellites. 

Each  specific  mission  where  cold  gas  thrusters  will 
be  used  places  a  unique  set  of  requirements  on  the 
system.  For  example,  micro-spacecraft  with  masses 
of  1  —  50  kg,  will  require  thrust  levels  for  typical 
orbital  maintenance  maneuvers  of  0.1  —  lOmN.  Mis¬ 
sions  which  require  precise  attitude  control,  where 
the  most  important  factor  is  a  very  low  thrust,  will 
require  thrust  levels  between  10  —  100/xN. 

To  obtain  these  low  thrust  values,  small  nozzle 
dimensions  and  low  chamber  pressures  are  usually 
used.  These  result,  in  the  throat  Reynolds  number, 
Ret  =  2771  / 7T Rt ,  on  the  order  of  10  to  103.  Here, 
the  quantity  m  is  the  mass  flow  rate,  Rt  is  the  throat 
radius,  and  jiq  is  the  propellant  gas  viscosity  at  the 
stagnation  chamber  temperature.  Due  to  the  low 
Reynolds  number,  viscous  losses  are  significant  in 
small  scale  nozzles. 

The  availability  of  micromachined  cold  gas 
thrusters  creates  new  possibilities  for  creating  effec¬ 
tive  low  thrust  propulsion  systems.  This  is  due  to 
a  possibility  of  manufacturing  nozzles  with  a  throat 
size  of  10pm  and  which  use  high  chamber  pressures 
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(1-10  atm).  This  permits  low  thrust  levels  for  rather 
high  Reynolds  numbers.  For  instance,  for  axisym- 
metric  nozzles,  with  a  fixed  area  ratio,  a  ten-fold 
change  in  the  linear  dimension  requires  a  100-fold 
increase  in  the  chamber  pressure  to  preserve  a  pre¬ 
scribed  thrust  level.  However,  the  throat  Reynolds 
number  increases  only  by  a  factor  of  ten,  and^the 
problem  of  taking  into  account  viscous  losses  re¬ 
mains. 

To  model  these  types  of  gas  flows,  Navier-Stokes 
(NS)  and  Direct  Simulation  Monte  Carlo  (DSMC) 
techniques  are  commonly  used.  The  weaknesses  of 
the  NS  approach  as  applied  to  rarefied  flows  are  well 
known  in  a  qualitative  sense.  At  sufficiently  low 
Ret ,  reliable  results  can  be  obtained  only  on  the  ba¬ 
sis  of  kinetic  approach,  for  example,  by  the  DSMC 
method1.  The  computational  cost  of  the  DSMC 
method,  however,  rapidly  increases  as  the  Reynolds 
number  increases. 

A  performance  evaluation  for  a  2D  micronozzle 
was  performed  using  the  NS  equations  in  Ref.  for 
Reynolds  numbers  Ret  =  500  —  10000  and  by  the 
DSMC3  method  for  Ret  =  1  -  160-  The  authors- 
used  a  combined  Navier-Stokes/DSMC  approach  to 
analyze  nozzle  flows  at  low  Reynolds  numbers.  To 
decrease  the  computational  cost,  the  DSMC  ^calcu¬ 
lation  was  performed  only  for  the  supersonic  sec¬ 
tion  of  the  nozzle.  The  development  of  numerical 
algorithms  of  the  DSMC  method  as  well  as  .  an  in¬ 
crease  of  computer  power  allows  the  calculation  of 
the  entire  nozzle  for  Reynolds  numbers  of  500  and 
higher3,  which  is  close  to  the  area  of  applicability  of 
the  Navier-Stokes  equations. 

The  objective  of  the  present  paper  is  to  compare 
local  fiowfield  and  integrated  performance  proper¬ 
ties  of  low-density  nozzles  exhausting  into  a  per¬ 
fect  vacuum  as  predicted  by  the  continuum  method 
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basedAa  finite  volume  approximation  of  NS  equa¬ 
tions,  and  the  DSMC  method  based  on  molecular 
gasdynamics.  The  calculations  were  performed  for  a 
conical  micronozzle  in  the  range  of  Reynolds  num¬ 
bers  Ret  =  130  -  1300  and  for  the  Rothe6  nozzle 
for  Ret  =  120  and  270  (corresponding  to  Rothe’s 
B  =  260  and  590).  These  comparisons  are  expected 
to  provide  design  guidance  for  nozzles  which  may  ex¬ 
perience  low  Reynolds  number  flow  due  to  low  stag¬ 
nation  pressure  or  due  to  small  characteristic  throat 
dimension  (micronozzles).  For  both  of  these  types 
of  nozzles,  viscous  and  rarefaction  effects  can  lead 
to  severe  performance  degradation  relative  to  ideal 
nozzle  predictions.  Since  the  NS  approach  is  typ¬ 
ically  much  more  computationally  efficient,  it  may 
be  preferable  to  DSMC  for  many  design  tasks  if  suf¬ 
ficient  confidence  in  its  ability  to  predict  the  perfor¬ 
mance  of  this  class  of  nozzles  can  be  established. 

Several  authors  have  previously  simulated  the 
Rothe  nozzle  using  both  simplified^  and  more 
complete®  forms  of  the  NS  equations  and  also 
DSMC,®  however,  these  focused  for  the  most  part 
on  comparison  with  local  measurements  made  on  the 
centerline  of  the  nozzle  diverging  section  rather  than 
integrated  performance  quantities  more  directly  use¬ 
ful  to  nozzle  design  studies.  It  will  be  shown  that 
the  NS  results  agree  rather  well  with  the  centerline 
data  even  at  relatively  low  operating  pressure  or  Retl 
however,  the  trends  in  specific  impulse  (or  thrust) 
begin  to  appeax  non-physical  at  these  same  pres¬ 
sures.  For  the  micronozzle,  no  experimental  data  are 
available,  and  only  comparisons  between  the  DSMC 
and  NS  results  are  made  as  a  function  of  operating 
pressure  and  nozzle  geometry.  The  lowest  pressure 
considered  leads  to  non-physical  behavior  of  the  NS- 
predicted  thrust.  Detailed  comparison  of  the  pre¬ 
dicted  flowfields  between  the  methods  facilitates  un¬ 
derstanding  of  the  deficiencies  in  the  NS  approach. 

2.  Nozzle  Geometry 


2.1.  Micronozzle 

A  conical  micronozzle  with  a  54.7°  or  30°  angle 
converging  section  and  15°  angle  diverging  section 
(Fig.  1)  used  for  the  computations.  The  nozzle  can 
be  routinely  fabricated  by  etching  silicon.  The  noz¬ 
zle  has  a  sharp  corner  at  the  throat,  an  exit  diameter 
of  81.3  \i m,  and  develops  a  thrust  of  1  mN  for  10  atm 
chamber  pressure  of  helium.  Chamber  and  wall  tem¬ 
peratures  were  To  =  Tw  =  297  K.  Chamber  pressure 
was  changed  from  10  to  1  atm  for  a  performance 
evaluation,  which  corresponded  to  a  change  in  the 
throat  Reynolds  number  from  1300  to  130. 


Fig.  1  Geometry  of  micronozzle 


2.2.  Rothe  nozzle 


Calculations  were  also  made  to  directly  compare 
with  the  experimental  data  of  Rothe.®  The  subsonic 
and  supersonic  portions  of  the  nozzle  are  cones  hav¬ 
ing  half-angles  30°  and  20°,  respectively,  with  lon¬ 
gitudinal  radii  of  curvature  at  the  throat  equal  to 
one-half  of  the  throat  radius.  The  test  gas  was  ni¬ 
trogen  at  T0  =  300  K.  The  throat  Reynolds  numbers 
were  Ret  =  270  and  120. 

3.  Numerical  approaches 


3.1.  Continuum  approach 


Continuum  solutions  are  generated  with  the  full 
Navier-Stokes  code  GASP*®.  The  nominal  grid  con¬ 
tains  257  axial  and  129  radial  points,  with  cluster¬ 
ing  in  both  directions  to  resolve  gradients  across  the 
throat  and  in  the  wall  boundary  layer.  The  inlet 
is  modeled  as  a  subsonic  Riemann  inflow  bound¬ 
ary  with  an  assumed  stagnation  number  density  and 
temperature,  and  Mach  number  based  on  the  ideal 
one-dimensional  nozzle  value  at  the  inlet  area  ratio. 
Fixed  stagnation  pressure  and  temperature  bound¬ 
ary  conditions  were  also  attempted,  but  were  aban¬ 
doned  due  to  code  instabilities^  First  order  extrap o- 
Qiation  is  used  at  the  outflow  boundary,  ^fiesehound- 
ary  conditions  are  usually  used  fSrHUZzle  flow  mod¬ 
eling.  It  is  assumed  that  the  use  of  extrapolation 
for  subsonic  portion  of  the  boundary  layer  at  the 
nozzle  exit  does  not  affect  nozzle  performance®.  For 
the  Rothe  nozzle  cases,  calculations  were  made  for 
both  adiabatic  and  constant  temperature  (Tw  =  To) 
and  no-slip  wall  boundary  conditions.  The  exper¬ 
imental  results  suggest  that  the  wall  is  adiabatic, 
however,  the  constant  wall  temperature  calculations 
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were  made  to  evaluate  the  effect  of  the  wall  bound¬ 
ary  conditions  on  the  flow  structure  and  integrated 
performance  properties.  For  the  micronozzle,  a  no¬ 
slip,  constant  temperature  (Tw= 297  K)  wall  was  as¬ 
sumed. 

Results  were  post-processed  to  extract  integrated 
performance,  local  flowfield,  and  wall  properties. 
The  axial  mass  flow  and  inviscid  momentum  flux 
(thrust)  were  interpolated  onto  each  axial  cell-face 
and  summed  at  each  axial  station  to  monitor  so¬ 
lution  quality  and  also  provide  a  local  measure  of 
nozzle  performance.  The  nominal  grid  shows  a  max¬ 
imum  axial  variation  of  mass  flow  of  below  one  per¬ 
cent  (Fig.  2).  Oscillations  in  total  mass  flow  are 
noticeable  at  the  throat,  due  presumably  to  the  in¬ 
terpolation  procedure,  and  at  each  boundary,  due  to 
interpolation  and  also  potentially  the  boundary  con¬ 
dition  treatment.  Momentum  flux  oscillations  are 
much  smaller.  The  specific  impulse  profiles  show 
negligible  oscillation.  Grid  sequencing  studies  con¬ 
firm  that  the  integrated  profiles  are  essentially  inde¬ 
pendent  of  grid  size  for  a  grid  of  129x32  points  or 
larger  (Table.  1). 

m/mj 


Fig.  2  Variation  of  axial  mass  flow  for  different 
nominal  grids  (micronozzle  at  10  atm  chamber  pres¬ 
sure) 


grid  size 

m/mi 

thrust,  mN 

Isp/Ispide 

257x129 

0.9313 

0.8746 

0.9840 

257x65 

0.9313 

0.8747 

0.9840 

257x33 

0.9316 

0.8752 

0.9840 

129x33 

0.9304 

0.8727 

0.9840 

65x33 

0.9631 

0.8992 

09835 

Table  1  The  influence  of  grid  size  on  the  integrated 
parameters  (micronozzle  at  10  atm  chamber  pres¬ 
sure) 


cells.  Each  background  cell  has  its  own  interaction 
parameter.  This  parameter  governs  the  local  col¬ 
lision  resolution.  The  adaptation  of  the  resolution 
to  local  gradients  and  local  mean  free  paths  occurs 
during  the  modeling  process.  Such  a  combined  us¬ 
age  of  cell  and  free  cell  schemes  made  it  possible  to 
achieve  adequate  spatial  resolution  in  the  entire  flow 
field  (see  Ref.  ^).  In  conjunction  with  the  correct 
collision  number  for  even  small  number  of  particles 
in  a  cell,  this  allows  one  to  simulate  the  flow  at  very 
small  Knudsen  numbers. 

2)  Adaptive  grids.  Two  different  grids  are  utilized 
in  the  code:  first  one  is  for  collision  procedure,  and 
second  is  for  gas  dynamic  parameters  sampling.  The 
utilization  of  two  different  grids  for  collisions  and 
macroparameters  is  rather  important  since  differ¬ 
ent  adaptation  algorithms  are  in  use  for  these  grids. 
Adaptation  of  collision  grids  aims  at  an  accordance 
of  cell  size  to  the  value  of  local  mean  free  path. 
The  macroparameter  grid  catches  generally  the  flow 
gradients  and  requires  normally  the  number  of  cells 
much  less  than  for  collisions. 

3) .  Spatial  weighting  factors .  The  most  severe 
practical  problem  associated  with  DSMC  calcula¬ 
tions  for  axisymmetric  flows  is  a  small  fraction  of 
the  molecules  that  are  located  near  the  axis.  This 
has  led  to  the  introduction  of  radial  weighting  fac¬ 
tors  such  that  a  molecule  located  far  from  the  axis 
represents  more  real  molecules  than  one  near  the 
axis. 


The  DSMC  code  SMILE  was  used  in  the  present 
work,  which  was  effectively  used  previously  for  jet 
flows  ^  and  nozzle  fiows^.  The  following  features  of 
the  code  are  important  for  computations  were  per¬ 
formed: 

1).  The  method  that  combines  majorant  cell  and 
free  cell  schemes ^  of  the  DSMC.  The  computational 
domain  is  divided  into  uniform  so-called  background 


4) .  Adaptive  decomposition  onto  subdomains  with 
different  time  steps.  The  strong  differences  in  the 
number  densities  and  flow  gradients  which  are  ob¬ 
served  near  to  the  throat,  and  those  observed  in  the 
near  vicinity  of  nozzle  exit  plane,  make  us  split  the 
computational  domain  into  subdomains  with  differ¬ 
ent  time  steps  (At).  This  splitting  is  performed  in 
the  course  of  computation  and  permits  to  compute 
the  nozzle  flow  faster.  An  approach  iWA*  = 
const  is  used  for  coupling  of  these  subdomains. 

5) .  Particle  doubling.  One  of  the  main  reasons  of 


3 


large  computational  cost  is  the  long  time  required  to 
establish  the  steady  flow.  The  cost  of  unsteady  stage 
exceeds  often  fifty  percent  of  the  whole  computation. 
To  overcome  this  problem,  a  special  procedure  of 
doubling  of  particles  is  used.  Using  the  doubling 
procedure  allows  one  to  compute  the  unsteady  stage 
rather  fast  at  the  successive  increasing  of  the  number 
of  particles. 

The  diffuse  law  with  complete  accommodation 
was  used  for  gas/surface  interaction,  and  the  VHS 
model^  for  intermolecular  collisions.  The  Larsen- 
Borgnakke  model^  with  temperature-dependent  ro¬ 
tational  collison  number  was  used  for  a  energy  ex¬ 
change  between  translational  and  rotational  molecu¬ 
lar  modes.  A  total  number  of  model  particles  in  the 
computations  was  1  *  106  -  3  •  106  (depending  on  the 
value  of  the  throat  Reynolds  number)  in  the  entire 
computational  domain  ABCD  (Fig.  1). 

Two  types  of  inlet  radial  velocity  profiles  (uni¬ 
form  and  Poiseuille  profiles)  were  considered  in  the 
present  DSMC  calculations.  For  the  nozzle  with  a 
54.7°  convergence  section  the  inlet  plane  was  located 
at  the  entrance  of  convergence  section  (lines  EF  in 
Fig.  1)  or  at  the  entrance  of  cylindrical  chanel  (line 
AD).  For  the  nozzle  with  a  30°  convergence  section, 
line  AD,  which  is  the  beginning  of  convergence  sec¬ 
tion,  was  also  used  as  an  inlet  plane.  The  particles 
enter  the  flow  domain  in  accordance  to  Maxwellian 
distribution  function  with  density,  temperature,  and 
velocity  obtained  from  isentropic  relations.  The 
mean  velocity  for  Poiseuille  profile  equaled  the  uni¬ 
form  flow  velocity. 

As  an  example,  Figure  3  shows  the  profiles  of  ax¬ 
ial  velocity  near  the  nozzle  throat  for  different  inlet 
plane  conditions  listed  in  Table  2.  It  is  seen  that  for 
all  cases  considered  the  mass  flow  rate  and  the  vac¬ 
uum  specific  impulse  coincide.  It  should  be  noted 
that  the  formulation  of  stable  and  physically  real¬ 
istic  inflow  conditions  at  a  subsonic  boundary  for 
both  NS  and  DSMC  techniques  is  a  current  research 
problem. 

4.  Results  and  Discussion 

4.1.  Micronozzle 

Calculations  of  micronozzle  flows  are  made  for 
the  different  chamber  pressures.  Figure  4  shows 
Mach  number  contours  obtained  by  two  methods 
for  po  =  10  atm.  It  is  seen  that  the  results  are  in 
good  agreement  in  the  throat  area  up  to  x  =  25^m 
(M  =  2.75).  Further  downstream  the  difference 
increases.  The  inflections  in  the  pressure  contours 
indicate  a  weak  compression  wave  forms  near  the 
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Fig.  3  Axial  velocity  profiles  for  different  inlet  con¬ 
ditions  (micronozzle  at  1  atm  chamber  pressure,  x  = 
10/im) 


throat  and  is  reflected  at  the  nozzle  symmetry  line 
(Fig.  5).  In  solving  the  Navier-Stokes  equations,  the 
outflow  boundary  of  the  computational  domain  coin¬ 
cided  with  the  nozzle  exit  plane,  and  the  parameters 
at  the  nozzle  exit  plane  were  determined  by  extrap¬ 
olation  from  inside.  These  boundary  conditions  re¬ 
sult  in  the  isoline  M  =  1  being  propagated  directly 
to  the  exit  plane.  In  the  kinetic  model,  a  strong 
pressure  gradient  induced  by  gas  exhaust  into  vac¬ 
uum  leads  to  flow  acceleration  near  the  nozzle  lip. 
and  the  isoline  M  —  1  is  terminated  on  the  nozzle 
lip.  The  difference  in  the  boundary  conditions  leads 
to  a  decrease  of  the  subsonic  section  of  the  bound¬ 
ary  layer  in  DSMC  simulations  and  to  a  decrease  of 
the  incident  angle  of  the  compression  wave.  Obvi¬ 
ously,  not  only  the  thickness  of  the  subsonic  section 
of  the  boundary  layer  decreases,  but  also  the  thick¬ 
ness  of  the  boundary  layer  itself,  which  leads  to  a 
greater  expansion  of  the  flow:  a  greater  Mach  num¬ 
ber  and  a  lower  pressure  (Fig.  5).  Thus,  even  for 
a  rather  small  thickness  of  the  subsonic  section  of 
the  boundary  layer  the  use  of  extrapolation  leads  to 
a  significant  difference  between  the  continuum  and 
kinetic  results.  These  local  features  are  known  to 
greatly  influence  the  near-exit  plume  or  contamina¬ 
tion  backflow  properties  of  the  nozzle. 

An  additional  factor  that  affects  the  flow  struc¬ 
ture  is  a  rather  large  value  of  the  slip  velocity  (~  400 
m/sec  at  the  nozzle  lip).  As  long  as  the  slip  veloc¬ 
ity  does  not  exceed  100  m/sec,  the  pressures  on  the 
wall  of  the  supersonic  section  of  the  nozzle  are  in 
fairly  good  agreement  (Fig.  6).  As  the  slip  velocity 
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case 

1 

2 

3 

4 

angle  of  convergence  section 

54.7 

54.7 

54.7 

30 

location  of  inlet  plane 

line  EF 

line  AD 

line  AD 

line  AD 

velocity  profile  at  inlet 

Poiseuille 

Poiseuille 

uniform 

uniform 

m,  kg /sec 

4.8421e-8 

4.8676e-8 

4.8394e-8 

4.8388e-8 

thrust,  mN 

0.073704 

0.074517 

0.074145 

0.074283 

7 sp ,  sec  _ 

155.22 

156.11 

156.23 

156.54 

,  Table  2  The  influence  of  inlet  conditions  on  the  micronozzle  performance  (1  atm  chamber  pressure) 
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Fig.  4  Mach  number  contours  (dashed  lines  - 
GASP,  solid  -  SMILE,  10  atm  chamber  pressure) 


Fig.  5  Pressure  isolines  (dashed  lines  -  GASP,  solid 
-  SMILE,  10  atm  chamber  pressure) 


increases  downstream,  the  pressures  differ  by  a  fac¬ 
tor  of  two.  Obviously,  the  effect  of  the  velocity  slip 
is  significant  only  inside  the  boundary  layer.  This 
is  clearly  shown  in  Fig.  7,  which  shows  the  profiles 
of  the  axial  and  radial  velocities  at  the  nozzle  exit 
plane.  In  the  core  flow  (R/Re  <  0.6)  both  numerical 
approaches  give  roughly  equal  velocities,  but  a  sub¬ 
stantial  difference  is  observed  inside  the  boundary 
layer. 

Two  criteria  are  usually  used  to  evaluate  the  area 


of  applicability  of  the  continuum  approach:  break¬ 
down  parameter  B  =  Iff  I  and  local  Knud- 


sen  number  Kn  =  p  |  ff | •  Af  is  the  Mach  number, 
p  is  the  density,  A  is  the  local  mean  free  path,  and 
s  is  the  distance  along  the  streamline.  For  the  case 


P/PO  Us,  m/sec 
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Fig.  6  Pressure  distribution  and  slip  velocity  along 
the  nozzle  wall  (chamber  pressure  is  10  atm) 


under  consideration  (po  =  10  atm),  the  breakdown 
parameter  and  local  Knudsen  number  contours  are 
shown  in  Figs.  8  and  9.  It  is  seen  that  the  critical 
value  B  =  0.05,  which  gives  the  border  of  the  area 
of  applicability  of  continuum  approach,  is  reached 
only  near  the  nozzle  lip.  It  can  be  concluded,  there¬ 
fore,  that  the  differences  observed  in  the  fields  of 
gasdynamic  parameters  obtained  using  different  ap¬ 
proaches  axe  explained  only  by  the  presence  of  slip  at 
the  wall  and  incorrect  outflow  boundary  conditions 
used  in  solving  the  Navier-Stokes  equations. 

The  situation  becomes  even  worse  when  these 
boundary  conditions  are  used  for  a  lower  chamber- 
pressure  (po  —  1  atm).  In  this  case,  the  subsonic 
section  of  the  boundary  layer  obtained  using  the 
Navier-Stokes  equations  is  more  than  30%  of  the 
exit  radius,  which  makes  a  two-fold  difference  with 
the  DSMC  result  (Fig.  10).  This  leads  to  the  fact 
that  the  Mach  number  in  the  nozzle  exit  does  not  ex¬ 
ceed  M  =  1.8,  whereas  the  kinetic  approach  yields 
significantly  higher  Mach  numbers  up  to  M  =  2.7. 


o 


Fig.  8  Breakdown  parameter  (10  atm  chamber  Fig.  II  Pressure  isolines  (dashed  lines  -  GASP, 
pressure,  SMILE)  solid  -  SMILE,  1  atm  chamber  pressure) 


For  this  micronozzle,  the  pressure  fields  obtained 
by  the  two  approaches  differ  rather  significantly,  ex¬ 
cept  in  the  vicinity  of  the  throat  (Fig.  11).  Obvi¬ 
ously,  the  continuum  approach  cannot  be  used  for 
such  a  low  pressure.  The  breakdown  parameter 
reaches  the  critical  value  B  =  0.05  approximately 
in  the  middle  of  the  nozzle  (Fig.  12).  Figure  13 
shows  the  nonequilibrium  that  arises  between  trans¬ 
lational  temperatures  based  on  the  thermal  compo¬ 
nent  of  motion  in  the  axial  and  radial  directions.  It 
is  seen  that  the  difference  is  observed  not  only  near 
the  nozzle  lip,  but  also  along  the  wall  of  the  diver¬ 
gence  section  of  the  nozzle. 


Fig.  9  Local  Knudsen  numbers  (10  atm  chamber 
pressure,  SMILE) 


Fig.  12  Breakdown  parameter  (SMILE,  1  atm 
chamber  pressure) 


Fig.  13  Nonequilibrium  of  temperature  Tz/Ty 
(SMILE,  1  atm  chamber  pressure) 
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Velocity  profiles  presented  in  Fig.  14  show  that 
the  boundary  layers  practically  fill  the  entire  noz¬ 
zle  exit,  and  the  presence  of  a  significant  slip  veloc¬ 
ity  (up  to  1000  m/sec)  leads  to  a  substantial  radial 
spreading  of  the  flow.  The  effect  of  a  high  pressure 
gradient  at  the  nozzle  exit  plane  extends  upstream 
to  the  throat  (Fig.  15).  This  demonstrates  once 
again  that  the  outflow  boundary  conditions  used  in 
the  Navier-Stokes  solution  are  invalid.  Additional 
Navier-Stokes  calculations  are  planned  to  address 
any  influence  that  inclusion  of  the  geometry  and 
flowfield  downstream  of  the  nozzle  exit  plane  has 
on  these  quantities. 


Fig.  14  Axial  U  and  radial  V  velocity  profiles  at 
nozzle  exit  plane  (1  atm  chamber  pressure) 


Note  that  the  use  of  the  density  distribution  along 
the  nozzle  axis  as  a  criterion  of  applicability  of  the 
continuum  approach  is  not  enough.  The  density  dis¬ 
tributions  along  the  nozzle  axis  for  the  examined  val¬ 
ues  of  chamber  pressure  (po  =  10  and  1  atm), 'which 
were  obtained  by  the  two  different  approaches  are 
in  good  agreement  as  shown  in  Fig.  16.  Neverthe¬ 
less,  the  Mach  number  distributions  along  the  nozzle 
axis  (Fig.  17)  testify  that  the  continuum  approach 
is  applicable  only  for  higher  chamber  pressures. 

The  nozzle  performance  calculated  by  the  DSMC 
method  is  shown  in  Fig.  18  for  a  wide  range  of 
Reynolds  numbers  Ret .  A  ten-fold  decrease  of  Ret 
decreases  the  coefficient  of  discharge  from  91%  to 
85%,  and  the  specific  impulse  from  95%  to  90%  of 
the  value  of  the  ideal  specific  impulse  obtained  from 
isentropic  relations.  The  calculations  performed  by 
GASP  yield  slightly  higher  values  of  the  coefficient 
of  discharge  and  specific  impulse  for  Ret  =  1300, 


Fig.  15  Pressure  distribution  and  slip  velocity 
along  nozzle  wall  (chamber  pressure  is  1  atm))  . 


P'Po 


Fig.  16  Density  distribution  along  nozzle  axis 


whereas  for  Ret  =  130  they  predict  a  lower  value  of 
mass  flow  rate.  The  NS  calculation  for  Ret  =  130 
gives  a  non-physical  ISp/Ispideai  >  1-  Fig*  19  shows 
the  change  in  Isp  along  the  nozzle  axis.  It  is  of  inter¬ 
est  to  note  that  Isp/Ispidial  >  1  is  observed  near  the 
nozzle  throat.  As  shown  in  fig.  19,  under  sufficiently 
low  Reynolds  number  operating  conditions,  the  large 
expansion  ratios  (long  diverging  section  lengths)  mo¬ 
tivated  by  ideal  flow  analyses  lead  to  excessive  vis¬ 
cous  losses  and  thus  a  peak  in  nozzle  thrust  occurs 
upstream  of  the  nozzle  exit.  The  potential'  reduc¬ 
tion  in  thruster  size  and  weight  indicate  one  obvi¬ 
ous  additional  advantage  of  a  shorter  nozzle.  Since 
the  flow  through  a  shorter  nozzle  will  exhibit  a  thin¬ 
ner  near-exit  plane  boundary  layer,  thruster-induced 
contamination  backflow  may  also  be  decreased. 
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Fig.  17  Mach  number  distribution  along  nozzle 
axis 


ro^ideal  Isp/Isp  ideal 


Fig.  18  Performance  evaluation  of  the  micronozzle 


4.2.  Rothe  nozzle 


The  experimental  data  obtained  in  Ref.  ®  were 
used  in  the  papers^’®’®  to  estimate  the  reliability 
of  results  of  numerical  simulation  of  low-Reynolds- 
number  nozzle  flows.  A  detailed  comparison  of  gas 
dynamic  characteristics  obtained  using  the  contin¬ 
uum  and  kinetic  approaches  for  the  Rothe  nozzle  was 
performed  in  the  present  work.  In  addition,  a  per¬ 
formance  evaluation  was  conducted  for  two  throat 
Reynolds  numbers. 

The  density  fields  for  Ret  =  270  are  presented  in 
Fig.  20.  Figures  21  and  22  show  the  density  distri¬ 
bution  along  the  nozzle  axis  and  in  a  cross-section 
near  the  nozzle  exit  plane.  It  is  seen  that  both  nu¬ 
merical  approaches  are  in  good  agreement  with  each 
other  and  with  experimental  values  of  density. 


isp 


Fig.  19  Specific  impulse  along  nozzle  axis 


p/po 


Fig.  21  Density  distribution  along  the  nozzle  axis 


P^Paxis 


Fig.  22  Density  profiles  neax  the  nozzle  exit 
(Xf  Rt  =  18.7,  Ret  =  270) 

A  comparison  of  temperature  fields  is  shown  in 
Fig.  23.  In  the  core  flow  both  approaches  predict 
close  value  of  the  translational  temperature.  A  dif¬ 
ference  is  observed  only  near  the  nozzle  wall  and 
in  the  vicinity  of  the  nozzle  lip.  A  comparison  of 
the  measured  value  of  the  rotational  temperature 
along  the  nozzle  axis  (Fig.  24)  with  numerical  pre¬ 
dictions  shows  their  good  agreement.  However,  the 
continuum  approach  gives  significantly  higher  values 
of  temperature  near  the  nozzle  wall  (Fig.  25). 


Fig.  23  Temperature  contours  (Rot he  nozzle, 
Ret  =  270,  dashed  lines  -  GASP,  solid  -  SMILE) 

The  differences  in  Mach  number  contours  are  simi¬ 
lar  to  those  in  temperature  flowfields  (Fig.  26).  They 
are  apparently  connected  with  the  position  of  the 
sonic  line.  The  outflow  extrapolation  bondary  condi¬ 
tions  extend  upstream  to  ~  10%.  The  performance 
characteristics  presented  in  Table  3  show  that  for 
Ret  =  270  the  difference  in  two  approaches  does  not 
exceed  2%. 

It  should  be  noted  that  despite  a  rather  low 
Reynolds  number  Ret  —  270  the  flow  characteris¬ 
tics  obtained  using  two  approaches  are  fairly  close. 
In  the  case  of  the  micronozzle,  even  for  Ret  =  1300 


Fig.  24  Temperature  distribution  along  the  nozzle 
axis  (Ret  —  270) 

rn0 


R/Rw 

Fig.  25  Temperature  profiles  near  nozzle  exit 
(X/Rt  =  18.7,  Ret  =  270,  adiabatic  wall) 

the  two  numerical  approaches  give  slightly  greater 
differences.  This  is  probably  related  to  the  nozzle 
geometry.  The  micronozzle  (Fig.  1)  has  a  compara¬ 
tively  short  divergence  section  (~  3.7  nozzle  throat 
diameters),  and  the  length  of  this  section  for  the 
Rothe  nozzle  equals  10.  Thus,  the  effect  of  outflow 
boundary  conditions  is  significantly  smaller  for  the 
Rothe  nozzle. 

It  was  already  noted  that  as  Ret  decreases  the  ef¬ 
fect  of  outflow  conditiohs  increases.  The  Mach  num¬ 
ber  field  shown  in  Fig.  27  demonstrates  that  for  the 
continuum  approach  the  boundary  layers  have 
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Fig.  26  Mach  number  contours  (Rothe  nozzle, 
Ret  =  270,  dashed  lines  -  GASP,  solid  -  SMILE) 


Fig.  27  Mach  number  contours  (Rothe  nozzle, 
Ret  =  120,  dashed  lines  -  GASP,  solid  -  SMILE) 


already  merged.  This  leads  to  a  decrease  of  the  Mach 
number  along  the  nozzle  axis  (closed  Mach  number 
contour).  In  the  kinetic  approach,  the  Mach  number 
increases  almost  monotonically  along  the  nozzle  axis. 

Since  in  the  experimental  paper  ®  the  test  gas  was 
nitrogen,  at  low  Reynolds  numbers  it  is  necessary  to 
take  into  account  different  relaxation  rates  of  trans¬ 
lational  and  rotational  energies.  In  the  continuum 
approach,  it  is  assumed  in  the  present  calculations 
that  the  translational  Ttrans  and  rotational  Trot  tem¬ 
peratures  are  equal.  The  rotational  nonequilibrium 
conturs  (1  -Trot/TtTans)  presented  in  Fig.  28  shows  a 
significant  difference  of  rotational  and  translational 
temperatures  near  the  nozzle  wall.  The  values  of 
'Rtrans  and  Trot  along  the  nozzle  axis,  which  were 
obtained  using  the  kinetic  approach,  are  in  good 
agreement  with  experimental  results,  whereas  the 
temperature  predicted  by  the  continuum  approach 
is  different  (Fig.  29).  This  difference  increases  down¬ 
stream  and  reaches  30%  near  the  nozzle  exit  plane. 

These  features  of  the  flow  associated  with  flow  rar¬ 
efaction  lead  to  a  significant  difference  in  the  pre¬ 
dicted  value  of  the  vacuum  specific  impulse  (see  Ta¬ 
ble  3).  It  reaches  6%,  which  is  three  times  greater 


Fig.  28  Rotational  nonequilibrium  (1  Trot/Ttrans , 
Ret  =  120,  SMILE) 


TT0 


X/Rt 


Fig.  29  Temperature  distribution  along  the  nozzle 
axis  (Ret  =  120) 


than  in  the  case  of  Rothe  nozzle  flow  for  Ret  =  270. 

5.  Conclusions 

Low-Reynolds-number  flows  have  been  studied 
for  two  types  of  nozzles  in  a  wide  range  of  throat 
Reynolds  numbers.  It  is  shown  that  the  use  of  ex¬ 
trapolation  conditions  as  the  boundary  conditions 
for  the  Navier-Stokes  equations  at  the  nozzle  exit 
plane  leads  to  overpredicting  of  the  vacuum  specific 
impulse  even  for  rather  high  Ret.  As  Ret  decreases, 
the  thickness  of  the  subsonic  portion  of  the  boundary' 
layer  increases,  and  the  use  of  extrapolation  condi¬ 
tions  can  lead  to  an  increase  of  the  vacuum  specific 
impulse  instead  of  its  decrease  because  of  viscous 
losses. 

A  computational  domain  extended  downstream  is 
used  in  the  DSMC  method,  and  vacuum  boundary 
conditions  are  set  on  it.  This  permits  a  correct  sim¬ 
ulation  of  the  flow  inside  the  nozzle  operating  in  the 
vacuum  of  the  space  environment  and  the  back  flow 
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/-Code 


GASP 

SMILE 

GASP 

SMILE 

GASP 

SMILE 

GASP 

J3MILE 


Ret 


270 

270 

270 

270 

120 

.120 

120 

120 


wall  condition 


adiabatic 
adiabatic 
Tw  =  300 AT 
Tw  =  300AT 
adiabatic 
adiabatic 
Tw  =  30.0jRT 
Tw  =  300 AT 


m,  kg/ sec 

thrust,  N 

Ap,  sec 

77l/ TTlideal 

1 $p  / 1 Spiral 

1.995e-05 

1.310e-02 

0.8981 

2.013e-05 

1.298e-02 

65.77 

0.9062 

0.8505 

1.983e-05 

1.331e-02 

68.50 

0.8927 

0.8858 

2.014e-05 

1.305e-02 

66.08 

0.9067 

0.8545 

8.364e-06 

5.324e-03 

64.96 

0.8539 

0.8400 

8.549e-06 

5.120e-03 

61.07 

0.8728 

0.7897 

8.276e-06 

5.489e-03 

67.68 

0.8450 

0.8752 

8.550e-06 

5.124e-03 

61.39 

0.8729 

0.7938 

Table  3  Performance  of  Rothe  nozzle 


expanding  around  the  nozzle  lip. 

In  the  case  of  micronozzles  with  a  comparatively 
short  supersonic  section  (~  3.7),  the  effect  of  incor¬ 
rect  outflow  extrapolation  conditions  is  quite  signifi¬ 
cant  and  is  observed  even  for  rather  high  Ret .  There¬ 
fore,  the  boundary  conditions  that  take  into  account 
a  subsonic  character  of  boundary  layer  flow  should 
be  used  in  the  Navier-Stokes  equations  to  evaluate 
the  performance  of  low-Reynolds-number  nozzles. 

Obviously,  additional  experimental  data  would 
be  useful  in  further  quantifying  the  accuracy  of 
these  numerical  methods.  To  this  end,  a  program 
has  been  initiated  to  measure  the  performance  of 
low  Reynolds  number  nozzle  and  micronozzle  flows. 
Based  on  available  instrumentation,  variations  in 
thrust  on  the  order  of  100/iN  and  in  specific  impulse 
on  the  order  of  0.5s  can  currently  be  resolved. 
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